Unrecognized diversity, genetic structuring, and phylogeography of the genus Triplophysa (Cypriniformes: Nemacheilidae) sheds light on two opposite colonization routes during Quaternary glaciation that occurred in the Qilian Mountains

Abstract In recent years, the influence of historical geological and climatic events on the evolution of flora and fauna in the Tibetan Plateau has been a hot research topic. The Qilian Mountain region is one of the most important sources of biodiversity on the Qinghai‐Tibet Plateau. Many species existed in the region during the Pleistocene glacial oscillation, and the complex geographical environment provided suitable conditions for the survival of local species. The shrinkage, expansion, and transfer of the distribution range and population size of species have significant effects on genetic diversity and intraspecific differentiation. To reveal the effects of geological uplift and climate oscillation on the evolution of fish populations in the Qilian Mountains, we investigated the genetic structure, phylogenetic relationship, and phylogeographical characteristics of genus Triplophysa species in the Qilian Mountains using the mitochondrial DNA gene (COI), three nuclear genes (RAG1, sRH, and Myh6) and 11 pairs of nuclear microsatellite markers. We collected 11 species of genus Triplophysa living in the Qilian Mountains, among which Triplophysa hsutschouensis and Triplophysa papillosolabiata are widely distributed in the rivers on the northern slope of the Qilian Mountains. There was a high degree of lineage differentiation among species, and the genetic diversity of endemic species was low. The different geographical groups of T. papillosolabiata presented some allogeneic adaptation and differentiation, which was closely related to the changes in the river system. Except for the population expansion event of T. hsutschouensis during the last glacial period of the Qinghai‐Tibet Plateau (0.025 MYA), the population sizes of other plateau loach species remained stable without significant population expansion. Starting from the east and west sides of the Qilian Mountains, T. hsutschouensis, and T. papillosolabiata showed two species colonization routes in opposite directions. The geological events of the uplift of the Qinghai‐Tibet Plateau and the climatic oscillation of the Quaternary glaciation had a great influence on the genetic structure of the plateau loach in the Qilian Mountains, which promoted the genetic differentiation of the plateau loach and formed some unique new species. The results of this study have important guiding significance for fish habitat protection in the Qilian Mountains.


| INTRODUC TI ON
Speciation is the result of a complex combination of factors that may accelerate or hinder speciation (Nosil et al., 2005). For example, different selection effects of different environmental conditions on populations can accelerate speciation (Wan et al., 2018;Wu et al., 2020;Zeng et al., 2020), and secondary contact of postglacial population dispersal can delay species differentiation (Hinojosa et al., 2019;Laurence, 2019). Geologic historical processes have a significant impact on species diversity distribution patterns (He et al., 2019;Zhang et al., 2019). Geomorphic changes caused by mountain uplift (Elias, 2009;Zhou et al., 2014) and climate change (Vila et al., 2005) will lead to population isolation, hinder genetic material exchange and increase differentiation among populations, which is also an important mechanism of speciation.
Biodiversity was strongly influenced by global climate and environmental changes during the Quaternary glaciation, with repeated climatic oscillations leading to the extinction of many phylogenetic lineages, followed by recolonization during interglacial periods (Hewitt, 2000;Taberlet et al., 1998). Some cold-tolerant species have been able to persist in glacial epoch and survive in refuges (Schmitt, 2007;Stewart et al., 2010), resulting in a high degree of species diversity in major refuge areas, whereas populations established during recolonization usually have low species diversity (Jablonski et al., 2016).
The unique geological features of the Tibetan Plateau make many areas a natural refuge for many species (Chen et al., 2017Wang et al., 2018;Zhou et al., 2014). Qilian Mountain, located at the northeastern edge of the Qinghai-Tibet Plateau, is a large and extremely important mountain range (Liu et al., 2002) that was uplifted along with the uplift of the Qinghai-Tibet Plateau (Li, 1963;Li & Fang, 1998;Li et al., 2001). Three inland rivers were formed in the northern area of the Qilian Mountains: the Shiyang, Hei, and Shule Rivers. In the southern region of the Qilian Mountains are the Datong and Huangshui Rivers, which are important tributaries of the Yellow River basin. There is some debate about the uplift time of the Qilian Mountains, with some scholars suggesting that the uplift time was between 10 and 20 million years ago (MYA) (George et al., 2001;Wang, 1997;Yue et al., 2001), while some scholars believe that the uplift time began approximately 8.3 MYA (Fang et al., 2005;Li, 1963). The geological uplift process of the Qilian Mountains experienced many drastic climatic changes and at least two large glacial periods (Liu, 1946(Liu, , 1951. These events had a great impact on the distribution pattern of species diversity from the Qilian Mountains (Zhao et al., 2011;Zhou et al., 2013). The biogeography of Schizopygopsis chilianensis (Li & Chang, 1974) (Cypriniformes: Cyprinidae), an endemic fish from the Qilian Mountains, suggests that this is a pattern of population diffusion along the Shiyang River system in the east of the Qilian Mountains to the west. In this process, species allopatronal evolution was related to the two great ice ages in the Qilian Mountains (Zhao et al., 2011).
Triplophysa Rendahl 1933 (Plateau loach) is a fish genus belonging to the family Nemacheilidae, included within the order Cypriniformes. This genus is widely distributed in the Tibetan Plateau and surrounding water systems and is a special group of loaches adapted to the cold environment of the Tibetan Plateau (Chen et al., 1996). Nine species of Plateau loach were first recorded in the Qilian Mountains (Zhu, 1989). With the discovery of new species and the redefinition of species, 11 species of plateau loach are reported in the Qilian Mountains (Wang, 1991;Wang, Yang, et al., 2015;Wu & Wu, 1992).
Due to the narrow distribution range and difficulties in obtaining samples, some Triplophysa species were rarely investigated for a long time in the Qilian Mountains (Feng et al., 2017;He et al., 2011;Wang et al., 2020). Are there any hidden species taxa neglected in the Qilian Mountains? Meanwhile, what about the distribution ranges and population genetic structure of these plateau loaches in this mountain system. These issues are still unknown .
In this study, we collected and analyzed data from the plateau loach fish (genus Triplophysa) populations in the Qilian Mountains, in order to: (1) estimate species diversity and species geographic distribution of the plateau loach fish; (2) reveal the implicit diversity of widely distributed species in the region; and (3) reconstruct the biogeographic history of the plateau loach fish population in the Qilian Mountains. Finally, we tested whether the species differentiation of the plateau loach in the Qilian Mountains was driven by specific geological events and whether the species distribution was likely to coincide with the changes in water systems on the genetic structure of the plateau loach in the Qilian Mountains, which promoted the genetic differentiation of the plateau loach and formed some unique new species.
The results of this study have important guiding significance for fish habitat protection in the Qilian Mountains.

K E Y W O R D S
Continental River, cryptic species, DNA barcode, genetic differentiation, plateau loach

T A X O N O M Y C L A S S I F I C A T I O N
Biodiversity ecology, Biogeography, Conservation genetics in the Qilian Mountains. Therefore, we expected to suppose that the genetic diversity of the Plateau loach was related to the geological changes in the region.

| Samplecollection
In this study, a total of 548 plateau loach fish samples were collected from the Datong, Huangshui, Shiyang, Hei, and Shule Rivers in the Qilian Mountains from the years of 2017 to 2020 ( Figure 1 and Table 1). As outgroup, nine Hedinichthys yarkandensis (Day 1877) were collected from the Shule River for further analysis. These specimens were captured using a cage net. Morphological identification of fishes is mainly based on taxonomic books (Kottelat, 2012;Wang, 1991;Zhu, 1989). The right pectoral fin strip of the specimen was cut and stored in 95% ethanol for total DNA extraction. The voucher specimen was fixed in 10% formaldehyde solution and stored in the fish exhibition room of Gansu Fisheries Research Institute.
MISA software (http://pgrc.ipk-gater sleben.de/misa/) was used to search microsatellite loci from the genomic data of T. siluroides (NCBI taxonomy ID: 422203). The search conditions were as follows: the length of the microsatellite repeat unit was 3-6 bp, the number of repeat units was >5, and the flanking sequence of microsatellite loci was >300 bp. Using genomic DNA from 5 individuals of each species as a template, 100 pairs of microsatellite primers were screened.
The fluorescence labeling technique was used to detect microsatellite labeling polymorphism, ABI 3730XL DNA analyzer (Thermo Fisher Scientific, USA) with a GeneScan-500Liz size standard was used for capillary electrophoresis, liz-500 was used as reference, GeneMarker V2.2.0 software (http://www.softg eneti cs.com) was used for microsatellite labeling typing. Eleven pairs of highly specific and polymorphic microsatellite primers were obtained for subsequent experimental analysis (Table 3). The PCR amplification system was the same as that of the previous COI gene fragment.  (Rozas et al., 2003) based on the mitochondrial COI gene sequence.

F I G U R E 1 Sample collection of the
Using Network 4.6 software (Polzin & Daneschmand, 2003), a Media-Joining haplotype Network map was constructed based on COI, RAG1, sRH1, and Myh6 gene sequences, and the connections between haplotypes were determined using maximum parsimony.
Using Hedinichthys yarkandensis as an outgroup, the phylogenetic tree was constructed by maximum likelihood (ML), Bayesian (BI), and neighbor-joining (NJ) methods. jModelTest Version 0.1.1 software (Posada, 2008) was used to select the optimal base substitution model for each partition. ML analysis was constructed using the RAxML 8.2.10 software package (Stamatakis, 2014 (Ronquist & Huelsenbeck, 2003) was used for Bayesian analysis. The starting tree was a random tree, and the optimal model was GTR + G. A total of 1.5 × 10 7 generations were run, and samples were taken every 1000 generations.
The first 3000 sampling trees (Burn-in = 3000) were discarded, the remaining trees were reserved to build a majority principle consensus tree, and the Bayesian Posterior Probability (BPP) of the evolutionary tree was calculated.
To analyze phylogeny in a time-calibrated framework, BEAST V.
2.4.6 software was used to estimate the relative time to the most recent common ancestor (Bouckaert et al., 2014). To avoid the unrealistic assumption of strict molecular clocks, we used a loose molecular clock model for analysis (Drummond et al., 2006). Two runs (

| Populationhistorydynamicanalysis
To infer population history dynamics, nucleotide mismatch analysis (Mismatch) was estimated for each species in DnaSP 4.2 software (Rozas et al., 2003). Arlequin V3.5 software was used for AMOVA and estimation of Tajima's D and Fu's Fs values in the neutral evolution test to judge whether the population was balanced and stable and to calculate the genetic differentiation index F st value between each two populations (repetitions 1000) (Excoffier et al., 2007;Schneider & Excoffier, 1999).
In BEAST 2.4.6, a coalesce-based Bayesian Skyline plot (BSP) was used to reconstruct the historical population size dynamics of each species (except T. leptosoma) (Drummond & Rambaut, 2007).
JModelTest Version 0.1.1 software (Posada, 2008) was run twice for each species under the BIC standard, and the simulation under the MCMC method ran 50 million generations (5000, burn-in = 25%).
ESS was evaluated for the stability of posterior branching probability in Tracer 1.6 software (effective value >200), and the outputs of the two runs were combined in the LogCombiner 2.4.6 module. The BSP diagram was made in Tracer 1.6 software.

The ancestral distribution of Triplophysa hsutschouensis and
Triplophysa papillosolabiata were reconstructed by RASP v4.2 (Reconstruct Ancestral State in Phylogenies, http://mnh.scu.edu.cn/ soft/blog/RASP). The population was divided into three distribution areas according to water system: A: Shiyang River system; B: Heihe River system; C: Shule River system. The Bayesian final tree was TA B L E 2 Primers of mitochondrial DNA gene and nuclear gene sequences in this study.  An individual-based Bayesian clustering method was used to detect the optimal genetic structure of the population in STRUCTURE 2.3.4 (Hubisz et al., 2009). Twenty independent simulations for each K value ranging from 1 to 20 were run, and one million generations were run at a time, discarding the first 100,000 generations. STRUCTURE HARVESTER (Earl & Vonholdt, 2012), a web-based software program, was used to parse and summarize the STRUCTURE output data and infer the optimal true population.

| Speciesdistribution
A total of 548 specimens of plateau loach Triplophysa were collected in the Qilian Mountains, and 11 species were identified.
T. shiyangensis and T. robusta had the highest haplotype diversity (0.889 and 0.915, respectively), while T. hsutschouensis had the lowest haplotype diversity (0.167). The highest nucleotide diversity was found in T. shiyangensis (0.00636), and the lowest was found in T.
The 11

| Geneticdifferentiation
The optimal structure model separated the genetic variation into 14 clusters, and the other acceptable models were 4 and 11 clusters, respectively ( Figure 2). The species differentiation in the phylogenetic tree based on mitochondrial genes and nuclear genes was in good agreement with the optimal cluster shown in the structure diagram based on microsatellite data (Figures 3 and 4). When the K value was 14, T. papillosolabiata with high genetic diversity was further di- The fifth population (CSC5) was composed of individuals from Shiyang River, but this population was divided into the sympatric distribution of T. wuweiensis in the optimal structure model.  (Table 6) of the five populations of T. papillosolabiata showed that 82.06% of the molecular differences were within the populations, and 17.94% of the molecular differences were between the populations. The genetic differentiation among the populations was also extremely significant (F st = 0.17944; p < .001). The genetic differentiation index analysis based on nuclear genes between species revealed significant genetic differentiation between 11 species pairs, except that T.
leptosoma and T. papillosolabiata shared the same haplotype based on sRH1 and Myh6 gene sequences (Table 7).

| Phylogenyandhaplotypenetworkevolution
Based on the COI gene sequence, a phylogenetic tree with consistent topological structure was obtained by using maximum likelihood (ML), Bayesian inference (BI), and neighbor-joining (NJ) methods, but the support values of some nodes were different.
Therefore, only the NJ phylogenetic tree is given in this study ( Figure 4). The phylogenetic tree shows that 11 plateau loach species in the Qilian Mountains possess unique haplotypes. In the topological structure of the phylogenetic tree, a single clade with high support was formed according to morphological species classification. Two clades were clearly identified from the phylogenetic tree. One clade was composed of three species, that is, T.

| Populationhistorydynamicsandspecies differentiation time
The analysis of population dynamics provides an important perspective for the evolutionary history of loach species in the Qilian Mountains ( Figure 7 and Table 8

| Speciesdistributionandcolonization
There has been no detailed investigation of the species distribution of plateau loach in the Qilian Mountains, and only a few sporadic investigations have been carried out. Some new endemic species of plateau loach have been discovered successively (Li & Chang, 1974;Zhao, 1984;Zhu & Wu, 1981), but the precise range of species remains to be revealed. In this study, based on an extensive specimen collection, we present a detailed distribution of all 11 species of plateau loach in the Qilian Mountains. Triplophysa siluroides and T.
scleroptera are distributed in the Datong River and Huangshui River tributaries of the Yellow River in the southern Qilian Mountains.
Triplophysa robusta and T. stoliczkae are mainly distributed in the Huangshui River but also exist in the Gulang River, a tributary of the Shiyang River, which may be the remains of the Wuwei Basin after its isolation from the Yellow River, corresponding to their low genetic diversity (Feng, 1963). When the species are distributed discretely or mixed with discrete and continuous distribution, the marginal population is in a relatively isolated state, and it is easy to produce a small population completely separated from the main population.
Such populations can evolve into new species over time (Coyne & Orr, 2004). Triplophysa wuweiensis is only distributed in the Xida River in the west of Shiyang River, and is considered to be a new species formed in the Miocene due to the geographical isolation between the Wuwei Basin and the Yellow River system (Li & Chang, 1974).
Triplophysa shiyangensis is distributed only in the Shiyang River system and is considered a specialized species independently adapted to the water environment of the Shiyang River (Zhao & Wang, 1983).
Triplophysa hsutschouensis is widely distributed in the Shiyang, Heihe, and Shule Rivers, and is regarded as an effective species because its body length is 7.18-7.60 times of body height. Wang (1991)  can reduce the waste of body temperature in cold areas (Wang, Yang, et al., 2015;Wang et al., 2016;Yang et al., 2021). This may be due to scale degradation during the Pleistocene Ice Age, which was an adaptive evolution of T. hsutschouensis to a cold environment. Different landscape features and climate patterns can drive directional selection (Sobel et al., 2010), which may play an important role in the speciation process of T.
hsutschouensis. The process of population expansion is often accompanied by drastic population genetic changes. The genetic diversity of this marginal population decreased due to bottleneck effect or genetic drift (Gao & Gao, 2016;Hirsch et al., 2015;Yang et al., 2016). According to the genetic diversity of geographical populations in different water systems (Table 1) Schizopygopsis chilianensis (Zhao et al., 2011).
Triplophysa tenuis, T. strauchii, and T. papillosolabiata are mainly distributed in the Tarim Basin (Wu & Wu, 1992). The inland river system of the Hexi Corridor is the eastern edge of the species distribution area. T. tenuis was mainly distributed in the Shule and Heihe Rivers, while T. strauchii was only collected in the Linze section of the middle reaches of the Hei River but not in the Shule River. This phenomenon of a discontinuous distribution area may be due to the influence of the ice age, and the Hei population of T. strauchii could take refuge in the lower reaches of the Hei River at low altitudes so that the population could survive. However, the lower reaches of the Shule River had high altitudes and could not form a shelter, leading to the extinction of the Shule River population. This phenomenon also exists in the species distribution of Rana chensinensis (Zhou et al., 2012). Triplophysa papillosolabiata is widely distributed in all water systems of the Hexi Corridor, which may be related to its strong ecological adaptation. It can not only feed in the rapids of the Dang River at an elevation of 2300 m but also survive in trickle water and even survive in some small water bodies frozen to the bottom in winter (Zhao, 1984). According to geographical distribution and genetic diversity, the PIC of the two population of T. papillosolabiata  (Dufresnes et al., 2016;Gómez & Lunt, 2007). As the intersection area of the ancient Yellow and Tarim River systems, the Qilian Mountain area was influenced by the geological uplift of the Qinghai-Tibet Plateau many times and formed the existing water system. In the course of strong geological structure and river development, T. hsutschouensis colonized from east to west, while T. papillosolabiata colonized from west to east along the northern slope of the Qilian Mountains. The effects of ice age and geographical isolation promoted the genetic differentiation of plateau loach and even formed some endemic species (Feng et al., 2017;Lv et al., 2020;Wu et al., 2013). Thus, these factors laid the foundation for the current pattern of species distribution.

F I G U R E 5
Map of minimal evolutionary network of haplotypes in the Triplophysa fishes. (The area of circles is proportional to the haplotype frequencies, and mv1-mv14 are missing haplotypes. Lines linking haplotypes indicate the evolutionary paths among haplotypes, vertical bars on the linking lines represent mutation steps between haplotypes).

| Crypticdiversity
Combined with the degree of delicacy in species differentiation of plateau loach Triplophysa Wang et al., 2020), the existence of a highly differentiated lineage of plateau loach in the Qilian Mountains is consistent with the existing morphological taxonomy (Zhu, 1989). In the phylogenetic tree of mtDNA, 11 species have formed their own phylogenetic clades, and T. leptosoma and T. papillosolabiata are sister groups with great morphological differences (Wang, 1991;Zhao & Wang, 1983). Triplophysa leptosome mainly lives in environments with fine-flowing water, while T. papillosolabiata generally lives in environments with deep water and slowflowing water. From the perspective of the living environment, they have their own unique ecological niches. As seen from the structure F I G U R E 6 Network haplotype evolution based on three nuclear genes. (The area of circles is proportional to the haplotype frequencies, and mv1-mv24 are missing haplotypes. Lines linking haplotypes indicate the evolutionary paths among haplotypes, vertical bars on the linking lines represent mutation steps between haplotypes).
diagram based on microsatellite data, T. leptosoma was differentiated from the Changma or Danghe populations of T. papillosolabiata, which also corresponds to the geographical location of the specimen collection. The MJN map of the nuclear genes showed that T. leptosoma and T. papillosolabiata shared the same haplotype in sRH1 and Myh6, while the two species had their own unique haplotype in the nuclear gene RAG1, which indicated that the differentiation time of T. leptosoma and T. papillosolabiata was relatively young.
As seen from the structure diagram, when the K value was 11, T. papillosolabiata could be divided into different groups according to water system (except that the Suganhu population merged into the population of the Hei River system), but T. hsutschouensis does not show a similar phenomenon, which may be related to the strong exotic evolution ability of the plateau loach (Zhao, 1984). When the

| Speciesdifferentiationanddrainagehistory
The genetic differentiation and species origin of plateau fish are closely related to the uplift and climate change of the Tibetan Plateau (He et al., 2006). It is important to consider the differentia-  Nanshan into the Hexi Corridor and eventually merged into the ancient Hei River, which was consistent with the geological analysis results (Feng, 1963). The Suganhu population may have been formed by the Yanglong population in the upper reaches of the Taolai River, which was connected with the Harteng River due to river capture in the interglacial period. The Harteng River is the main water source of Sugan Lake. It is generally believed that there are three relatively independent river systems in the Hexi Corridor, namely, the Shiyang, Hei, and Shule River systems, and the Sugan Lake and Harteng River systems are merged into the Shule River system. According to the population differentiation of T. papillosolabiata, the Sugan Lake system was never historically connected with the Shule River system but was connected with the Taolai River in the upper reaches of the Heihe River.

| Implicationsforconservation
Population genetic diversity is an important part of biodiversity.
The level, formation mechanism, and distribution pattern of genetic diversity can reveal the origin and evolution of species and are also the basis of species conservation. The maintenance of genetic diversity prevents the loss of the evolutionary potential of species (Bartáková et al., 2019;). The genetic differentiation characteristics of the plateau loach fish in the Qilian Mountains play an important role not only in fish conservation but also in water ecological stability in the Qilian Mountains of China. According to the principle of setting evolutionarily significant units (ESUs) and management units (MUs) (Frankham, 2003;Moritz, 1994), our analysis shows that due to the unique genetic characteristics of T. shiyangensis and T.

| CON CLUS IONS
In this study, we evaluated the genetic differentiation of the genus

ACK N OWLED G M ENTS
We are grateful to Ka W, Zhang YS, and Hu YB for help to obtain the samples.

F I G U R E 11
Reconstructed ancestral distribution of Triplophysa hsutschouensis (I) and Triplophysa papillosolabiata (II) based on BBM. (a: Shiyang River system; b: Heihe River system; c: Shule River system; * represent other ancestral ranges. The different colors of the nodes represent different ancestral ranges. The number in the ring represents the node's coded. The values at the node represent support values in Bayesian analysis, and the posterior probability of less than 50 is not shown).